{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "from scipy.optimize import bisect"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f44f7150",
   "metadata": {},
   "source": [
    "# Model\n",
    "\n",
    "Preferences are given by : \n",
    "\n",
    "$$ U(c,h) = \\frac{c^{1-\\sigma_c}}{1-\\sigma_c} + \\phi \\frac{h^{1-\\sigma_h}}{1-\\sigma_h} $$\n",
    "\n",
    "The production function is \n",
    "\n",
    "$$h = \\alpha_1 + \\alpha_0 m$$\n",
    "\n",
    "The budget constraint is such that \n",
    "\n",
    "$$y = c + p m$$\n",
    "\n",
    "There are two regions, us and eu. Upon substituting $h$, we have the following problem: \n",
    "\n",
    "$$ \\max_{m} \\frac{(y-pm)^{1-\\sigma_c}}{1-\\sigma_c} + \\phi \\frac{(\\alpha_1 + \\alpha_0m)^{1-\\sigma_h}}{1-\\sigma_h} $$\n",
    "\n",
    "The FOC of this problem is \n",
    "\n",
    "$$ p(y - pm)^{-\\sigma_c} = \\alpha_0 \\phi (\\alpha_1 + \\alpha_0 m)^{-\\sigma_h} $$\n",
    "\n",
    "We first program up the solution, using bisection (no derivative involved)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "id": "01fd442f",
   "metadata": {},
   "outputs": [],
   "source": [
    "def foc(m,sigma_c,sigma_h,phi,alphas,y,p):\n",
    "\treturn p*(y-p*m)**(-sigma_c) - alphas[0]*phi*(alphas[1]+alphas[0]*m)**(-sigma_h)\n",
    "def optm(sigma_c,sigma_h,phi,alphas,y,p):\n",
    "\treturn bisect(foc,0.01*y,y*0.95,args=(sigma_c,sigma_h,phi,alphas,y,p))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "87cca3d1",
   "metadata": {},
   "source": [
    "These are the parameters that we end up getting, but we check here that the function works. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "id": "72da1031",
   "metadata": {},
   "outputs": [],
   "source": [
    "alphas_us = [6.66,0]\n",
    "sigma_c = 2\n",
    "sigma_h = 3\n",
    "phi = 0.207"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "98fafde7",
   "metadata": {},
   "source": [
    "We start with income and price normalized in U.S. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "12e2648f",
   "metadata": {},
   "outputs": [],
   "source": [
    "p_us = 1\n",
    "y_us = 1"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e2e61684",
   "metadata": {},
   "source": [
    "We can check that the function returns s in the U.S. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "id": "73f4bbf7",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.14995735503423932"
      ]
     },
     "execution_count": 15,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "optm(sigma_c,sigma_h,phi,alphas_us,y_us,p_us)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "id": "5bafd2b2",
   "metadata": {},
   "outputs": [],
   "source": [
    "s_us = 0.15\n",
    "h_us = 1.0"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "id": "97827789",
   "metadata": {},
   "outputs": [],
   "source": [
    "alphas_eu = [6.66,0.22]\n",
    "p_eu = 0.61\n",
    "y_eu = 0.78"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "id": "17f16d7b",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.09622527852617895"
      ]
     },
     "execution_count": 18,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "p_eu*optm(sigma_c,sigma_h,phi,alphas_eu,y_eu,p_eu)/y_eu"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "id": "b139b8ea",
   "metadata": {},
   "outputs": [],
   "source": [
    "s_eu = 0.096\n",
    "h_eu = 1.047"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "id": "376e26d9",
   "metadata": {},
   "outputs": [],
   "source": [
    "def uopt(sigma_c,sigma_h,phi,alphas,y,p):\n",
    "\tm = optm(sigma_c,sigma_h,phi,alphas,y,p)\n",
    "\treturn ((y-p*m)**(1-sigma_c))/(1-sigma_c) + phi*((alphas[1]+alphas[0]*m)**(1-sigma_h))/(1-sigma_h) \n",
    "def prod(alphas,m):\n",
    "\treturn alphas[1] + alphas[0]*m\n",
    "def hopt(sigma_c,sigma_h,phi,alphas,y,p):\n",
    "\tm = optm(sigma_c,sigma_h,phi,alphas,y,p)\n",
    "\treturn prod(alphas,m)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "id": "5691c041",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "1.0394607817832697"
      ]
     },
     "execution_count": 23,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "hopt(sigma_c,sigma_h,phi,alphas_eu,y_eu,p_eu)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e52353c8",
   "metadata": {},
   "source": [
    "# Calibration\n",
    "\n",
    "We start the exercise with finding U.S. parameters that fit the moments. We are solving a set of equations for two unknowns, $\\alpha_0$ and $\\phi$. We use the methods of moments to minimize the squared distance of moments. We use a derivative free algorithm "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "id": "a9f55003",
   "metadata": {},
   "outputs": [],
   "source": [
    "def criterion_us(theta,sigma_c,sigma_h,alphas,y,p,moms):\n",
    "\talphas[0] = theta[0]\n",
    "\tphi = theta[1]\n",
    "\tm = optm(sigma_c,sigma_h,phi,alphas,y,p)\n",
    "\tmoms_s = np.zeros(2)\n",
    "\tmoms_s[0] = p*m/y\n",
    "\tmoms_s[1] = prod(alphas,m)\n",
    "\treturn (moms[0] - moms_s[0])**2 + (moms[1] - moms_s[1])**2"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "209b7224",
   "metadata": {},
   "source": [
    "These are the numbers presented in Section 2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 29,
   "id": "f419d91b",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(6.666657902257493, 0.20761274012255118)"
      ]
     },
     "execution_count": 29,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "from scipy.optimize import minimize\n",
    "\n",
    "opt_us = minimize(criterion_us,[5.0,0.3],args=(sigma_c,sigma_h,[0.0,0.0],y_us,p_us,[s_us,h_us]))\n",
    "alphas_us[0] = opt_us.x[0]\n",
    "phi = opt_us.x[1]\n",
    "alphas_us[0],phi\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "5eacecdf",
   "metadata": {},
   "source": [
    "In Europe, we seek to estimate $p_{EU}$ and $\\alpha_{1,EU}$ using $(s_{EU},h_{EU})$ as moments"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 30,
   "id": "bcf770d6",
   "metadata": {},
   "outputs": [],
   "source": [
    "def criterion_eu(theta,sigma_c,sigma_h,phi,alphas,y,moms):\n",
    "\talphas[1] = theta[0]\n",
    "\tp = theta[1]\n",
    "\tm = optm(sigma_c,sigma_h,phi,alphas,y,p)\n",
    "\tmoms_s = np.zeros(2)\n",
    "\tmoms_s[0] = p*m/y\n",
    "\tmoms_s[1] = prod(alphas,m)\n",
    "\treturn (moms[0] - moms_s[0])**2 + (moms[1] - moms_s[1])**2"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "df417188",
   "metadata": {},
   "source": [
    "These are the numbers reported in Section 2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 33,
   "id": "4936d7b9",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(0.21439909063306478, 0.599575021285484)"
      ]
     },
     "execution_count": 33,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "alphas_eu[0] = alphas_us[0]\n",
    "opt_eu = minimize(criterion_eu,[0.4,0.8],args=(sigma_c,sigma_h,phi,alphas_eu,y_eu,[s_eu,h_eu]))\n",
    "alphas_eu[1] = opt_eu.x[0]\n",
    "p_eu = opt_eu.x[1]\n",
    "alphas_eu[1],p_eu"
   ]
  }
 ],
 "metadata": {
  "interpreter": {
   "hash": "cf2a50979671a58939829e6829efb726aa5da11149213b77bd50351f899d04fb"
  },
  "kernelspec": {
   "display_name": "Python 3.8.5 64-bit ('base': conda)",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.8.12"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}
